install.packages('tinytex')
library(tinytex)
tlmgr_search('framed.sty')  # 搜索包含 framed.sty 文件的 LaTeX 包
tlmgr_install('framed')     # 安装 framed 包
tlmgr_update()
tinytex::install_tinytex()
citation("sbgcop")
48844+30625
library(countrycode)
countrycode("PHL", "iso3c", "country.name")
countrycode("BRA", "iso3c", "country.name")
log(7/10)
log(10/7)
load("/Users/chongchen/Downloads/dataverse_files (8)/coefs_2fe.rda")
library(MCMCpack)
install.packages("MCMCpack")
library(MCMCpack)
m01 <- as.mcmc(coefs$coef_01)
m01
# R options
options(
htmltools.dir.version = FALSE, # for blogdown
show.signif.stars = FALSE,     # for regression output
warm = 1
)
# Set dpi and height for images
library(knitr)
opts_chunk$set(tidy = FALSE, fig.height = 3, dpi = 300, message=FALSE, warning=FALSE, cache=TRUE)
# ggplot2 color palette with gray
color_palette <- list(gray = "#999999",
salmon = "#E69F00",
lightblue = "#56B4E9",
green = "#009E73",
yellow = "#F0E442",
darkblue = "#0072B2",
red = "#D55E00",
purple = "#CC79A7")
# For nonsese...
#library(emo)
htmltools::tagList(rmarkdown::html_dependency_font_awesome())
library(ggplot2)
theme_set(theme_bw())
library(dlstats)
cran_stats(c("ggplot2", "sf", "tmap", "sp","spdep", "rgdal")) %>%
ggplot(aes(end, downloads, group=package, color=package)) +
geom_line() + geom_point(aes(shape=package)) + xlab("") +
scale_x_date(date_breaks = "6 month", date_labels =  "%Y-%m") + theme(text = element_text(size=8))
library(dplyr)
library(ggplot2)
library(dlstats)
cran_stats(c("ggplot2", "sf", "tmap", "sp","spdep", "rgdal")) %>%
ggplot(aes(end, downloads, group=package, color=package)) +
geom_line() + geom_point(aes(shape=package)) + xlab("") +
scale_x_date(date_breaks = "6 month", date_labels =  "%Y-%m") + theme(text = element_text(size=8))
setwd("~/Dropbox/Civil war alliances/rebelalliance/Analysis/ISQ_RR/Replication_ISQ")
rm(list=ls())
# Function to load packages
loadPkg=function(toLoad){
for(lib in toLoad){
if(! lib %in% installed.packages()[,1])
{install.packages(lib, repos='http://cran.rstudio.com/')}
suppressMessages( library(lib, character.only=TRUE))}}
# Load libraries
packs=c("readr", "readxl", "dplyr", "countrycode", "ggplot2", "amen","statnet","VGAM",
"amen", "purrr", "texreg", "abind", "ergm","reshape2",
'parallel','foreach','doParallel',"haven", "readstata13")
loadPkg(packs)
#### set working directory to the folder named "Replication_ISQ"
#setwd("Replication_ISQ")
#####################
source("Code/helpFUN.R")
source("Code/function_interaction_plot.R")
load("Data/glm_df4.RData")
glm_df <- glm_df4
rm(glm_df4)
## creade dyad anme
glm_df <- glm_df %>%
dplyr::mutate(dyad = paste(RBLside_A, RBLside_B, sep = "-"))
# create lag Y
glm_df <- glm_df %>%
dplyr::arrange(dyad, year) %>%
dplyr::group_by(dyad) %>%
dplyr::mutate(lag_y = dplyr::lag(alliance, n=1))
library(spduration)
#transform to spdur format
glm_df <- as.data.frame(glm_df)
glm_df <- add_duration(glm_df, "alliance", unitID = "dyad", tID ="year", freq = "year", ongoing = FALSE)
## create no-alliance years polynomial terms
glm_df <- glm_df  %>%
dplyr::mutate(t_yr = duration,
t_yrsq = duration^2,
t_yrcub = duration^3)
logit_f2 <- glm(alliance ~ co_ideology + allCOETH +splinter_indirect2+
rebels_count + pop_log + gdppc_log + milper + rugg_prop + post_cold +
t_yr + t_yrsq + t_yrcub,
data=glm_df, family='binomial'(link = "logit"))
logit_f3 <- glm(alliance ~ co_constituent +splinter_indirect2+
rebels_count + pop_log + gdppc_log + milper + rugg_prop + post_cold +
t_yr + t_yrsq + t_yrcub,
data=glm_df, family='binomial'(link = "logit"))
logit_f_coef <- scaled_coef_cluster(ModelResults = list(logit_f2), data = glm_df,
clusterid = "ccode", subvar = TRUE, vars = c(
"(Intercept)",
#"t_yr", "t_yrsq","t_yrcub",
"post_cold", "pop_log","gdppc_log",
"milper","rebels_count","rugg_prop", "splinter_indirect2",
"allCOETH","co_ideology"
#,"co_constituent"
))
logit_f_coef <- logit_f_coef +
theme(legend.position = "none") +scale_x_discrete(labels=c("Intercept",
#  "t",
# expression(t^2),
#  expression(t^3),
'Post-Cold War',
"Population",
'GDP per capita',
'Military personnel',
"Number of rebel groups",
"Terrain ruggedness",
"Splinter",
"Co-ethnic",
"Co-ideological"#,
# "Co-constituent"
)) +
ggtitle("")
library(postregplots)
logit_f2_mar1 <- density_ridgeplot(logit_f2, n.sim = 1000, varname = c('allCOETH'),
data = glm_df, val1 = 0, val2 = 1, clusterid = "ccode")
logit_f2_mar1 <-logit_f2_mar1 + xlab("Average of First Difference") +
scale_y_discrete(expand = c(-0.1, 0),labels=c("")) +
ggtitle("Co-ethnic") +
theme(axis.text = element_text(size=14),
text = element_text(size=14),
plot.title = element_text(size = 18))
mean(logit_f2_mar1$data$value) #0.07076735
## marginal for model2 and model 3
logit_f2_mar <- density_ridgeplot(logit_f2, n.sim = 1000, varname = c('co_ideology'),
data = glm_df, val1 = 0, val2 = 1, clusterid = "ccode")
mean(logit_f2_mar$data$value) #0.1615785
logit_f2_mar <- logit_f2_mar + xlab("Average of First Difference") +
ggtitle("Co-ideological") +
scale_y_discrete(expand = c(-0.1, 0),labels=c("")) +
theme(axis.text = element_text(size=14),text = element_text(size=14),
plot.title = element_text(size = 18))
library(grid)
library(gridExtra)
f <- grid.arrange(logit_f2_mar,logit_f2_mar1, ncol = 2)
f
f
grid.arrange(logit_f2_mar,logit_f2_mar1, ncol = 2)
## predicted
d = predict_probs(ModelResults = logit_f2, n.sim = 1000, varname = c('co_ideology'),
data = glm_df, val1 = 0, val2 = 1,intervals=1, clusterid = "ccode")
mean(logit_f2_mar1$data$value) #0.07076735
## for multinomial
data <- glm_df
table(data$dv_nominal)
data$dv_nominal <- factor(data$dv_nominal, labels = c("noalliance", "informal", "formal"))
table(data$dv_nominal)
library(VGAM)
nominal_tfit2 <- vglm(dv_nominal ~ co_ideology + allCOETH + splinter_indirect2+
rebels_count + pop_log + gdppc_log + milper + rugg_prop + post_cold +
t_yr + t_yrsq + t_yrcub,
data= data, multinomial(refLevel="noalliance"))
nominal_tfit3 <- vglm(dv_nominal ~ co_constituent + splinter_indirect2+
rebels_count + pop_log + gdppc_log + milper + rugg_prop + post_cold +
t_yr + t_yrsq + t_yrcub,
data= data, multinomial(refLevel="noalliance"))
coef_mlogit_t <- mlogit_scaled_coef(fit = list(nominal_tfit2), subvar = TRUE, vars = c(
"(Intercept)",
#"t_yr", "t_yrsq","t_yrcub",
"post_cold", "pop_log","gdppc_log",
"milper","rebels_count","rugg_prop","splinter_indirect2",
"allCOETH","co_ideology"
#, "co_constituent"
))
coef_mlogit_t <- coef_mlogit_t + theme(legend.position = "none")+ scale_x_discrete(labels=c("Intercept",
#  "t",
# expression(t^2),
#expression(t^3),
'Post-Cold War',
"Population",
'GDP per capita',
'Military personnel',
"Number of rebel groups",
"Terrain ruggedness",
"Splinter",
"Co-ethnic",
"Co-ideological"
)) +
ggtitle("")
p_ideology <- vglm_sim(mod = nominal_tfit2, var = "co_ideology", n_sim = 1000)
p_ideology$data %>% dplyr::group_by(type) %>% dplyr::summarise(mean(Pred))
p_ideology$data %>%group_by(type) %>% summarise(mean(Pred))
#1 formal        -0.0230
#2 informal       0.103
p3 <- p_allCOETH + scale_y_discrete(expand = c(-0.1, 0),
labels=c("Formal alliance","Informal alliance")) +
xlab("Average of First Difference") +
ggtitle("Co-ethnic") +
theme(axis.text = element_text(size=12),
text = element_text(size=12),
plot.title = element_text(hjust = .5, size = 14))
p_allCOETH <- vglm_sim(mod = nominal_tfit2, var = "allCOETH", n_sim = 1000)
p_allCOETH$data %>% dplyr::group_by(type) %>% dplyr::summarise(mean(Pred))
#1 formal        -0.0230
#2 informal       0.103
p3 <- p_allCOETH + scale_y_discrete(expand = c(-0.1, 0),
labels=c("Formal alliance","Informal alliance")) +
xlab("Average of First Difference") +
ggtitle("Co-ethnic") +
theme(axis.text = element_text(size=12),
text = element_text(size=12),
plot.title = element_text(hjust = .5, size = 14))
grid.arrange(p4, p3, nrow = 1)
p4 <- p_ideology + scale_y_discrete(expand = c(-0.1, 0),
labels=c("Formal alliance","Informal alliance")) +
xlab("Average of First Difference") +
ggtitle("Co-ideological") +
theme(axis.text = element_text(size=12),
text = element_text(size=12),
plot.title = element_text(hjust = .5, size = 14))
grid.arrange(p4, p3, nrow = 1)
rm(list=ls())
# Function to load packages
loadPkg=function(toLoad){
for(lib in toLoad){
if(! lib %in% installed.packages()[,1])
{install.packages(lib, repos='http://cran.rstudio.com/')}
suppressMessages( library(lib, character.only=TRUE))}}
# Load libraries
packs=c("readr", "readxl", "dplyr", "countrycode", "ggplot2", "amen","statnet","VGAM",
"amen", "purrr", "texreg", "abind", "ergm","reshape2",
'parallel','foreach','doParallel',"haven", "readstata13")
loadPkg(packs)
#### set working directory to the folder named "Replication_ISQ"
#setwd("Replication_ISQ")
#####################
source("Code/helpFUN.R")
source("Code/function_interaction_plot.R")
load("Data/glm_df4.RData")
glm_df <- glm_df4
rm(glm_df4)
## creade dyad anme
glm_df <- glm_df %>%
dplyr::mutate(dyad = paste(RBLside_A, RBLside_B, sep = "-"))
# create lag Y
glm_df <- glm_df %>%
dplyr::arrange(dyad, year) %>%
dplyr::group_by(dyad) %>%
dplyr::mutate(lag_y = dplyr::lag(alliance, n=1))
library(spduration)
#transform to spdur format
glm_df <- as.data.frame(glm_df)
glm_df <- add_duration(glm_df, "alliance", unitID = "dyad", tID ="year", freq = "year", ongoing = FALSE)
## create no-alliance years polynomial terms
glm_df <- glm_df  %>%
dplyr::mutate(t_yr = duration,
t_yrsq = duration^2,
t_yrcub = duration^3)
##############################################################################
# A17: Robustness Check: Logistic Regression Results for Alliance Using an Alternative Shared Ideological Constituency Variable (1946-2015)
##############################################################################
glm_df <- glm_df %>%
dplyr::mutate(co_ethnic_nuance1 = ifelse(ideolnationlist ==1 & allCOETH ==1, 1, 0),
co_ethnic_nuance2 = ifelse(allCOETH == 1 & (RBL_A_forge_ideolnat == 0 & RBL_B_forge_ideolnat == 0), 1, 0)) %>%
dplyr::mutate(co_ethnic_alternative = ifelse(co_ethnic_nuance1 == 1 | co_ethnic_nuance2 == 1, 1, 0))
table(glm_df$co_ethnic_nuance1)
table(glm_df$co_ethnic_nuance2)
table(glm_df$co_ethnic_alternative)
logit_f2_1 <- glm(alliance ~ co_ideology_new1 + allCOETH +splinter_indirect2+
rebels_count + pop_log + gdppc_log + milper + rugg_prop + post_cold +
t_yr + t_yrsq + t_yrcub,
data=glm_df, family='binomial'(link = "logit"))
logit_f_coef <- scaled_coef_cluster(ModelResults = list(logit_f2_1 ), data = glm_df,
clusterid = "ccode", subvar = TRUE, vars = c(
"(Intercept)",
# "t_yr", "t_yrsq","t_yrcub",
"post_cold", "pop_log","gdppc_log",
"milper","rebels_count","rugg_prop", "splinter_indirect2",
"allCOETH","co_ideology_new1"
))
logit_f_coef <- logit_f_coef +  theme(legend.position = "none")+  scale_x_discrete(labels=c("Intercept",
# "t",
#  expression(t^2),
#  expression(t^3),
'Post-Cold War',
"Population",
'GDP per capita',
'Military personnel',
"Number of rebel groups",
"Terrain ruggedness",
"Splinter",
"Co-ethnic",
"Co-ideological")) +
ggtitle("")
library(postregplots)
logit_f2_mar1 <- density_ridgeplot(logit_f2_1, n.sim = 1000, varname = c('allCOETH'),
data = glm_df, val1 = 0, val2 = 1, clusterid = "ccode")
logit_f2_mar1 <-logit_f2_mar1 + xlab("Average of First Difference") +
scale_y_discrete(expand = c(-0.1, 0),labels=c("")) +
ggtitle("Co-ethnic") +
theme(axis.text = element_text(size=14),
text = element_text(size=14),
plot.title = element_text(size = 18))
mean(logit_f2_mar1$data$value) #[1] 0.06407264
## marginal for model2 and model 3
logit_f2_mar <- density_ridgeplot(logit_f2_1, n.sim = 1000, varname = c('co_ideology_new1'),
data = glm_df, val1 = 0, val2 = 1, clusterid = "ccode")
logit_f2_mar <- logit_f2_mar + xlab("Average of First Difference") +
ggtitle("Co-ideological") +
scale_y_discrete(expand = c(-0.1, 0),labels=c("")) +
theme(axis.text = element_text(size=14),text = element_text(size=14),
plot.title = element_text(size = 18))
mean(logit_f2_mar$data$value) #[1] 0.1517364
grid.arrange(logit_f2_mar,logit_f2_mar1, ncol = 2)
##############################################################################
# A18: Robustness Check: Multinomial Logistic Regression Results for Formal and Informal Alliance Using an Alternative Shared Ideological Constituency Variable (1946-2015)
##############################################################################
data <- data %>%
dplyr::mutate(co_ethnic_nuance1 = ifelse(ideolnationlist ==1 & allCOETH ==1, 1, 0),
co_ethnic_nuance2 = ifelse(allCOETH == 1 & (RBL_A_forge_ideolnat == 0 & RBL_B_forge_ideolnat == 0), 1, 0)) %>%
dplyr::mutate(co_ethnic_alternative = ifelse(co_ethnic_nuance1 == 1 | co_ethnic_nuance2 == 1, 1, 0))
nominal_tfit2_1 <- vglm(dv_nominal ~ co_ideology_new1 + allCOETH + splinter_indirect2+
rebels_count + pop_log + gdppc_log + milper + rugg_prop + post_cold +
t_yr + t_yrsq + t_yrcub,
data= data, multinomial(refLevel="noalliance"))
coef_mlogit_t <- mlogit_scaled_coef(fit = list(nominal_tfit2_1), subvar = TRUE, vars = c(
"(Intercept)",
# "t_yr", "t_yrsq","t_yrcub",
"post_cold", "pop_log","gdppc_log",
"milper","rebels_count","rugg_prop","splinter_indirect2",
"allCOETH","co_ideology_new1"
))
coef_mlogit_t <- coef_mlogit_t + theme(legend.position = "none")+ scale_x_discrete(labels=c("Intercept",
#   "t",
#  expression(t^2),
#  expression(t^3),
'Post-Cold War',
"Population",
'GDP per capita',
'Military personnel',
"Number of rebel groups",
"Terrain ruggedness",
"Splinter",
"Co-ethnic",
"Co-ideological")) +
ggtitle("")
p_ideology <- vglm_sim(mod = nominal_tfit2_1, var = "co_ideology_new1", n_sim = 1000)
##############################################################################
# A18: Robustness Check: Multinomial Logistic Regression Results for Formal and Informal Alliance Using an Alternative Shared Ideological Constituency Variable (1946-2015)
##############################################################################
data <- data %>%
dplyr::mutate(co_ethnic_nuance1 = ifelse(ideolnationlist ==1 & allCOETH ==1, 1, 0),
co_ethnic_nuance2 = ifelse(allCOETH == 1 & (RBL_A_forge_ideolnat == 0 & RBL_B_forge_ideolnat == 0), 1, 0)) %>%
dplyr::mutate(co_ethnic_alternative = ifelse(co_ethnic_nuance1 == 1 | co_ethnic_nuance2 == 1, 1, 0))
## for multinomial
data <- glm_df
table(data$dv_nominal)
data$dv_nominal <- factor(data$dv_nominal, labels = c("noalliance", "informal", "formal"))
table(data$dv_nominal)
library(VGAM)
##############################################################################
# A18: Robustness Check: Multinomial Logistic Regression Results for Formal and Informal Alliance Using an Alternative Shared Ideological Constituency Variable (1946-2015)
##############################################################################
data <- data %>%
dplyr::mutate(co_ethnic_nuance1 = ifelse(ideolnationlist ==1 & allCOETH ==1, 1, 0),
co_ethnic_nuance2 = ifelse(allCOETH == 1 & (RBL_A_forge_ideolnat == 0 & RBL_B_forge_ideolnat == 0), 1, 0)) %>%
dplyr::mutate(co_ethnic_alternative = ifelse(co_ethnic_nuance1 == 1 | co_ethnic_nuance2 == 1, 1, 0))
nominal_tfit2_1 <- vglm(dv_nominal ~ co_ideology_new1 + allCOETH + splinter_indirect2+
rebels_count + pop_log + gdppc_log + milper + rugg_prop + post_cold +
t_yr + t_yrsq + t_yrcub,
data= data, multinomial(refLevel="noalliance"))
coef_mlogit_t <- mlogit_scaled_coef(fit = list(nominal_tfit2_1), subvar = TRUE, vars = c(
"(Intercept)",
# "t_yr", "t_yrsq","t_yrcub",
"post_cold", "pop_log","gdppc_log",
"milper","rebels_count","rugg_prop","splinter_indirect2",
"allCOETH","co_ideology_new1"
))
coef_mlogit_t <- coef_mlogit_t + theme(legend.position = "none")+ scale_x_discrete(labels=c("Intercept",
#   "t",
#  expression(t^2),
#  expression(t^3),
'Post-Cold War',
"Population",
'GDP per capita',
'Military personnel',
"Number of rebel groups",
"Terrain ruggedness",
"Splinter",
"Co-ethnic",
"Co-ideological")) +
ggtitle("")
p_ideology <- vglm_sim(mod = nominal_tfit2_1, var = "co_ideology_new1", n_sim = 1000)
p4 <- p_ideology + scale_y_discrete(expand = c(-0.1, 0),
labels=c("Formal alliance","Informal alliance")) +
xlab("Average of First Difference") +
ggtitle("Co-ideological") +
theme(axis.text = element_text(size=14),
text = element_text(size=14),
plot.title = element_text(hjust = .5, size = 16))
p_ideology$data %>% group_by(type) %>% summarise(mean = mean(Pred))
##
p_allCOETH <- vglm_sim(mod = nominal_tfit2_1, var = "allCOETH", n_sim = 1000)
p3 <- p_allCOETH + scale_y_discrete(expand = c(-0.1, 0),
labels=c("Formal alliance","Informal alliance")) +
xlab("Average of First Difference") +
ggtitle("Co-ethnic") +
theme(axis.text = element_text(size=14),
text = element_text(size=14),
plot.title = element_text(hjust = .5, size = 16))
p_allCOETH$data %>% group_by(type) %>% summarise(mean = mean(Pred))
p_ideology$data %>% dplyr::group_by(type) %>% dplyr::summarise(mean = mean(Pred))
p_allCOETH$data %>% dplyr::group_by(type) %>% dplyr::summarise(mean = mean(Pred))
grid.arrange(p4, p3, nrow = 1)
mean(logit_f2_mar$data$value) #[1] 0.1517364
list.dirs("code")
list.files("code")
p_ideology$data %>% group_by(type) %>% summarise(mean = mean(Pred))
p_allCOETH$data %>% group_by(type) %>% summarise(mean = mean(Pred))
mean(logit_f2_mar$data$value) #[1] 0.1517364
mean(logit_f2_mar$data$value) #[1] 0.1517364
mean(logit_f2_mar1$data$value) #[1] 0.06407264
mean(logit_f2_mar$data$value) #0.1615785
rm(list=ls())
# Function to load packages
loadPkg=function(toLoad){
for(lib in toLoad){
if(! lib %in% installed.packages()[,1])
{install.packages(lib, repos='http://cran.rstudio.com/')}
suppressMessages( library(lib, character.only=TRUE))}}
# Load libraries
packs=c("readr", "readxl", "dplyr", "countrycode", "ggplot2", "amen","statnet","VGAM",
"amen", "purrr", "texreg", "abind", "ergm","reshape2",
'parallel','foreach','doParallel',"haven", "readstata13")
loadPkg(packs)
### set working directory to the folder named "Replication_ISQ"
#setwd("Replication_ISQ")
#####################
source("Code/helpFUN.R")
source("Code/function_interaction_plot.R")
###############################################################################################
### This is the code for the models in the main text: Figure 3-5
###############################################################################################
load("Data/glm_df4.RData")
glm_df <- glm_df4
rm(glm_df4)
## creade dyad anme
glm_df <- glm_df %>%
dplyr::mutate(dyad = paste(RBLside_A, RBLside_B, sep = "-"))
# create lag Y
glm_df <- glm_df %>%
dplyr::arrange(dyad, year) %>%
dplyr::group_by(dyad) %>%
dplyr::mutate(lag_y = dplyr::lag(alliance, n=1))
library(spduration)
#transform to spdur format
glm_df <- as.data.frame(glm_df)
glm_df <- add_duration(glm_df, "alliance", unitID = "dyad", tID ="year", freq = "year", ongoing = FALSE)
## create no-alliance years polynomial terms
glm_df <- glm_df  %>%
dplyr::mutate(t_yr = duration,
t_yrsq = duration^2,
t_yrcub = duration^3)
logit_f2 <- glm(alliance ~ co_ideology + allCOETH +splinter_indirect2+
rebels_count + pop_log + gdppc_log + milper + rugg_prop + post_cold +
t_yr + t_yrsq + t_yrcub,
data=glm_df, family='binomial'(link = "logit"))
logit_f3 <- glm(alliance ~ co_constituent +splinter_indirect2+
rebels_count + pop_log + gdppc_log + milper + rugg_prop + post_cold +
t_yr + t_yrsq + t_yrcub,
data=glm_df, family='binomial'(link = "logit"))
logit_f_coef <- scaled_coef_cluster(ModelResults = list(logit_f2), data = glm_df,
clusterid = "ccode", subvar = TRUE, vars = c(
"(Intercept)",
#"t_yr", "t_yrsq","t_yrcub",
"post_cold", "pop_log","gdppc_log",
"milper","rebels_count","rugg_prop", "splinter_indirect2",
"allCOETH","co_ideology"
#,"co_constituent"
))
logit_f_coef <- logit_f_coef +
theme(legend.position = "none") +scale_x_discrete(labels=c("Intercept",
#  "t",
# expression(t^2),
#  expression(t^3),
'Post-Cold War',
"Population",
'GDP per capita',
'Military personnel',
"Number of rebel groups",
"Terrain ruggedness",
"Splinter",
"Co-ethnic",
"Co-ideological"#,
# "Co-constituent"
)) +
ggtitle("")
library(postregplots)
logit_f2_mar1 <- density_ridgeplot(logit_f2, n.sim = 1000, varname = c('allCOETH'),
data = glm_df, val1 = 0, val2 = 1, clusterid = "ccode")
logit_f2_mar1 <-logit_f2_mar1 + xlab("Average of First Difference") +
scale_y_discrete(expand = c(-0.1, 0),labels=c("")) +
ggtitle("Co-ethnic") +
theme(axis.text = element_text(size=14),
text = element_text(size=14),
plot.title = element_text(size = 18))
mean(logit_f2_mar1$data$value) #0.07076735
list.files("code")
52/12
